; docformat = 'rst'
;
;+
;
; :Purpose:
;   Create 3D field netCDF file. Vertical levels are for mid-points of sigma levels, 
;     such that for the ith level:
;     p[i]=hybridA[i] + hybridB[i]*PS
;   This routine also checks that the grid is orientated correctly.
;
; :Inputs:
;   lonInput: longitude
;   latInput: latitude
;   input: 3D data in dimensions [lon, lat, lev, time]
;   time: array of Julian dates corresponding to each element in the input array
;   hybridAInput: Hybrid A values, dimension=[lev], Units=Pa
;   hybridBInput: Hybrid B values, dimension=[lev], Unitless
;   PSInput: surface pressure at each time step, dimensions=[lon, lat, time], units=Pa
;   fileName: Name of netCDF file to be output
;   title: string detailing the contents of the file, used for the global "Title" attribute
;   fieldName: name for 3D field in netCDF file
;   fieldUnits: units of 3D field
;
; :Keywords:
;   author: string to be used for global "Author" attribute
;
; :Outputs:
;   Outputs 1 if successful
;
; :Example::
;   
;   IDL> success=gauge_3D_write_netcdf(lon, lat, time, hya, hyb, PS, input, $
;     "/my/directory/N2O_loss_photolysis_climatology.nc", $
;     "N2O J-values for stratospheric photolysis", $
;     "J(N2O)", $
;     "cm-3 s-1", $
;     author="Matt Rigby")
;   
; :History:
; 	Written by: Matt Rigby, University of Bristol, Jun 7, 2013
;
;-
function gauge_3D_write_netcdf, lonInput, latInput, timeInput, hybridAInput, hybridBInput, PSInput, fieldInput, $
    fileName, title, fieldName, fieldUnits, $
    author=author

  compile_opt idl2
;  on_error, 2

  time=timeInput
  lon=lonInput
  lat=latInput
  ps=psInput
  hybridA=hybridAInput
  hybridB=hybridBInput
  field=fieldInput

  ;Check array dimensions
  if n_elements(field[*, 0, 0, 0]) ne n_elements(lon) then begin
    print, "NETCDF WRITE ERROR: Check input field longitude dimension (should be 1st dimension)"
    return, 0
  endif
  if n_elements(field[0, *, 0, 0]) ne n_elements(lat) then begin
    print, "NETCDF WRITE ERROR: Check input field latitude dimension (should be 2nd dimension)"
    return, 0
  endif
  if n_elements(field[0, 0, *, 0]) ne n_elements(hybridA) or $
    n_elements(field[0, 0, *, 0]) ne n_elements(hybridB) then begin
    print, "NETCDF WRITE ERROR: Check input field vertical dimension (should be 3rd dimension)"
    return, 0
  endif
  if n_elements(field[0, 0, 0, *]) ne n_elements(time) then begin
    print, "NETCDF WRITE ERROR: Check input field time dimension (should be 4th dimension)"
    return, 0
  endif

  ;Check surface pressure dimensions
  if n_elements(ps[*, 0, 0]) ne n_elements(lon) then begin
    print, "NETCDF WRITE ERROR: Check urface pressure longitude dimension (should be 1st dimension)"
    return, 0
  endif
  if n_elements(ps[0, *, 0]) ne n_elements(lat) then begin
    print, "NETCDF WRITE ERROR: Check surface pressure latitude dimension (should be 2nd dimension)"
    return, 0
  endif
  if n_elements(ps[0, 0, *]) ne n_elements(time) then begin
    print, "NETCDF WRITE ERROR: Check surface pressure time dimension (should be 3th dimension)"
    return, 0
  endif


  ;Check that longitudes are ordered correctly
  whW=where(lon gt 180., countW, complement=whE)
  if countW gt 0 then begin
    print, "Longitudes should be ordered from -180 to 180, re-orientating"
    fieldSize=size(field)
    fieldNew=fltarr(fieldSize[1:fieldSize[0]])
    fieldNew[0:countW-1, *, *, *]=field[whW, *, *, *]
    fieldNew[countW:-1, *, *, *]=field[whE, *, *, *]
    field=fieldNew

    psSize=size(ps)
    psNew=fltarr(psSize[1:psSize[0]])
    psNew[0:countW-1, *, *]=ps[whW, *, *]
    psNew[countW:-1, *, *]=ps[whE, *, *]
    ps=psNew

    lon=[lon[whW] - 360., lon[whE]]
  endif

  ;Check that latitudes are ordered correctly
  if lat[1] lt lat[0] then begin
    print, "Latitudes should be ordered from -90 to 90, re-orientating"
    lat=reverse(lat)
    ps=reverse(ps, 2)
    field=reverse(field, 2)
  endif

  ;Check that levels are ordered from TOA to surface
  pTest=hybridA + 100000.*hybridB
  if pTest[-1] lt pTest[0] then begin
    print, 'Re-orientating vertical levels (top of atmosphere to surface)"
    hybridA=reverse(hybridA)
    hybridB=reverse(hybridB)
    field=reverse(field, 3)
    pTest=reverse(pTest)
  endif
  print, 'Vertical levels from approximately: ' + string(pTest[0]/100., format='(f8.2)') + $
    ' to ' + string(pTest[-1]/100., format='(f8.2)') + 'hPa'
    
  ;Convert time
  timeType=size(time, /type)
  case 1 of    
    ;Julian dates
    (timeType eq 4 or timeType eq 5): begin
      caldat, time, month, day, year, hour, minute, second
      date=long(year)*10000L + long(month)*100L + long(day)
      dateSec=long(hour)*60L*60L + long(minute)*60L + long(second)
    end
    ;Integer: assume this is a month in a monthly climatology
    (timeType eq 2 or timeType eq 3): begin
      date=long(time)
      datesec=replicate(0L, n_elements(time))
    end
  endcase


  ;Write netCDF file
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  fID=ncdf_create(fileName, /clobber)
  
  ; Create required dimensions
  lonDID  = NCDF_DIMDEF(fID, 'lon', long(n_elements(lon)))
  latDID  = NCDF_DIMDEF(fID, 'lat', long(n_elements(lat)))
  levDID  = NCDF_DIMDEF(fID, 'lev', long(n_elements(HybridA)))
  timeDID  = NCDF_DIMDEF(fID, 'time', /unlimited)
  
  ; Create a variable to hold the data
  lonID = NCDF_VARDEF(fId,'lon', [LonDID], /float)
  latID = NCDF_VARDEF(fId,'lat', [LatDID], /float)
  hyamID = NCDF_VARDEF(fId,'hyam', [levDID], /float)
  hybmID = NCDF_VARDEF(fId,'hybm', [levDID], /float)
  psID= NCDF_VARDEF(fid, 'PS', [lonDID, latDID, timeDID], /float)
  dateID = NCDF_VARDEF(fId,'date', [TimeDID], /LONG)
  dateSecID = NCDF_VARDEF(fId,'datesec', [TimeDID], /LONG)
  fieldID = NCDF_VARDEF(fId, fieldName, [LonDID, LatDID, levDID, TimeDID], /float)
  
  ; Create some attributes
  NCDF_ATTPUT, fid, /Global, "Title", title
  if keyword_set(author) then begin
    NCDF_ATTPUT, fid, /Global, "Author", author
  endif
  
  NCDF_ATTPUT, fId, lonId, "long_name", "longitude"
  NCDF_ATTPUT, fId, lonId, "units", "degrees_east"
  
  NCDF_ATTPUT, fId, latId, "long_name", "latitude"
  NCDF_ATTPUT, fId, latId, "units", "degrees_north"

  NCDF_ATTPUT, fId, hyamId, "long_name", "Hybrid A coordinate"
  NCDF_ATTPUT, fId, hyamId, "units", "Pa"

  NCDF_ATTPUT, fId, hyamId, "long_name", "Hybrid B coordinate"
  NCDF_ATTPUT, fId, hyamId, "units", "None"

  NCDF_ATTPUT, fId, psID, "long_name", "Surface pressure"
  NCDF_ATTPUT, fId, psId, "units", "Pa"
  
  NCDF_ATTPUT, fId, fieldId, "long_name", fieldName
  NCDF_ATTPUT, fId, fieldId, "units", fieldUnits

  NCDF_ATTPUT, fId, dateID, "long_name", "Date"
  if date[0] gt 10000000L then begin
    NCDF_ATTPUT, fId, dateID, "units", "YYYYMMDD"    
  endif else begin
    NCDF_ATTPUT, fId, dateID, "units", "MM"    
  endelse

  NCDF_ATTPUT, fId, dateSecID, "long_name", "Seconds past 0000 UTC"
  NCDF_ATTPUT, fId, dateSecID, "units", "seconds"
  
  ; Leave definition mode and enter data write mode
  NCDF_CONTROL, fId, /ENDEF
  
  ; Write the data
  NCDF_VARPUT, fId, lonId, lon
  NCDF_VARPUT, fId, latId, lat
  NCDF_VARPUT, fId, hyamId, hybridA
  NCDF_VARPUT, fId, hybmId, hybridB
  NCDF_VARPUT, fId, psId, ps  
  NCDF_VARPUT, fId, dateID, date
  NCDF_VARPUT, fId, datesecID, datesec
  NCDF_VARPUT, fId, fieldId, field
  
  ; Done
  NCDF_CLOSE, fId

  return, 1

end